home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1997 August / Walnut Creek CDROM.7z / VOL_400 / 446_01 / DOC / FDMSTART / EX8 / CONSLAW.C < prev    next >
Encoding:
C/C++ Source or Header  |  1996-04-18  |  6.8 KB  |  275 lines
  1.                                                                  // -*- C++ -*-
  2. // $Id: $
  3. #include <ConsLaw.h>
  4. #include <MenuSystem.h>
  5. #include <createHCLSchemes.h>
  6. #include <createFluxFunc.h>
  7. #include <createInitFunc.h>
  8. #include <DrawFD.h>
  9.  
  10. ConsLaw:: ConsLaw ()
  11. {
  12.   TRACE("ConsLaw::ConsLaw")
  13.  
  14.   p_value = 0;  // gives p=1 field
  15. }
  16.  
  17.  
  18. BooLean ConsLaw:: ok () const
  19. {
  20.   TRACE("ConsLaw::ok")
  21.  
  22.   if (!u.ok())
  23.     errorFP("ConsLaw::ok","u is not rebinded to a field (empty handle)");
  24.   else if (!u->ok())
  25.       errorFP("ConsLaw::ok","u field is not ok");
  26.  
  27.   if (!u_prev.ok())
  28.     errorFP("ConsLaw::ok","u_prev is not rebinded to a field (empty handle)");
  29.   else if (!p->ok())
  30.       errorFP("ConsLaw::ok","u_prev field is not ok");
  31.  
  32.   if (!scheme.ok())
  33.     errorFP("ConsLaw::ok","scheme is not rebinded to a field (empty handle)");
  34.  
  35.   if (!flux.ok())
  36.     errorFP("ConsLaw::ok","flux is not rebinded to a function (empty handle)");
  37.  
  38.   if (!p.ok())
  39.     errorFP("ConsLaw::ok","p is not rebinded to a field (empty handle)");
  40.   else if (!p->ok())
  41.       errorFP("ConsLaw::ok","p field is not ok");
  42.  
  43.   if (time_points_for_plot.getNoMembers() < 1)
  44.     warningFP("ConsLaw::ok","no plots will be made! (no time points given)");
  45.  
  46.   return dpTRUE;
  47. }
  48.  
  49. void ConsLaw:: adm (MenuSystem& menu)
  50. {
  51.   TRACE("ConsLaw::adm")
  52.  
  53.   define (global_menu);   // define menu items
  54.   global_menu.prompt();   // prompt the user (read input data)
  55.   scan (global_menu);     // load user's input data into the class variables
  56. }
  57.  
  58. void ConsLaw:: solveProblem ()
  59. {
  60.   TRACE("ConsLaw::driver")
  61.  
  62.   if (!ok())
  63.     fatalerrorFP("ConsLaw::driver",
  64.          "object is not properly initialized, have you called adm?");
  65.   setIC ();            // set initial conditions
  66.   timeLoop ();            // perform time stepping
  67. }
  68.  
  69. void ConsLaw:: define (MenuSystem& menu, int level)
  70. {
  71.   menu.addItem (level,
  72.                 "grid",
  73.                 "grid",
  74.                 "GridLattice::scan syntax",
  75.                 "d=1, domain: [0,1] index: [0:20]",
  76.                 "S",
  77.                 'g', 'g');
  78.   menu.addItem (level,
  79.                 "time parameters",
  80.                 "timeprm",
  81.                 "TimePrm::scan syntax",
  82.                 "dt=0.05, t in [0,0.9]",
  83.                 "S",
  84.                 't','t');
  85.   menu.addItem (level,
  86.                 "time points for plot",
  87.                 "plot_times",
  88.                 "Set::scan syntax, eoi: <CR>",
  89.                 "0.2 0.5 0.8",
  90.                 "S", // improvement: could be intelligent from t in [0:?]
  91.                 'p','p');
  92.   menu.addItem (level,
  93.                 "scheme",
  94.                 "scheme",
  95.                 "classname in HCLSchemes hierarchy",
  96.                 "UpwindScheme1D",
  97.                 validationString(hierHCLSchemes()),
  98.                 's','s');
  99.   menu.addItem (level,
  100.                 "flux function",
  101.                 "flux",
  102.                 "functor, name in FluxFunc hierarchy",
  103.                 "FluxLinear",
  104.                 validationString(hierFluxFunc()),
  105.                 'f','f');
  106.   menu.addItem (level,
  107.                 "initial function",
  108.                 "u0",
  109.                 "functor, name in InitFunc hierarchy",
  110.                 "InitConst1",
  111.                 validationString(hierInitFunc()),
  112.                 '0','0');
  113.   menu.addItem (level,
  114.                 "uL",
  115.                 "uL",
  116.                 "value of u(x,t) at the left boundary",
  117.                 "1",
  118.                 "R1",
  119.                 'l','l');
  120.   menu.addItem (level,
  121.         "p-mean",
  122.         "E[p]",
  123.         "mean value of the p-field (f(u;p))",
  124.         "0",
  125.         "R1",
  126.         'P','P');
  127. }
  128.  
  129. void ConsLaw:: scan (MenuSystem& menu)
  130. {
  131.   TRACE("ConsLaw::scan")
  132.  
  133.   grid.rebind (new GridLattice(1));
  134.   grid->scan (menu.get("grid"));
  135.  
  136.   u.rebind (new FieldFD (grid(), "u"));
  137.   u_prev.rebind (new FieldFD (grid(), "u_prev"));
  138.  
  139.   tip.scan (menu.get("time parameters"));
  140.  
  141.   // add the end mark ';' to the end (required by TimePrm::scan):
  142.   time_points_for_plot.scan (menu.get("time points for plot") + " ;");
  143.  
  144.   // flux term:
  145.   flux_tp = menu.get("flux function");
  146.   flux.rebind (createFluxFunc(flux_tp,*this));
  147.  
  148.   // initial function
  149.   u0_tp = menu.get("initial function");
  150.   u0.rebind (createInitFunc(u0_tp,*this));
  151.  
  152.   // x=xMin boundary condition:
  153.   uL = menu.get("uL").getReal();
  154.   u->valueIndex(0) = u_prev->valueIndex(0) = uL;
  155.  
  156.   // scheme type:
  157.   scheme_tp = menu.get("scheme");
  158.   scheme.rebind (createHCLSchemes (scheme_tp, *this));
  159.  
  160.   // p-field, set p = exp(p_value):
  161.   p.rebind (new FieldFD (grid(),"p"));
  162.   p_value = menu.get("p-mean").getReal();
  163.   p->fill (exp(p_value));
  164.  
  165.   // plot file name:
  166.   file.open (CaseName);
  167.  
  168.   // write input data for a check:
  169. }
  170.  
  171. void ConsLaw:: setIC ()
  172. {
  173.   TRACE("ConsLaw::setIC")
  174.  
  175.   u_prev() = u0();
  176.   u_prev->valueIndex(grid->getBase(1)) = uL;
  177.  
  178.   // --- find the suitable dt ---
  179.   real dt;
  180.   real pmin, pmax;
  181.   p->minmax (pmin, pmax);
  182.  
  183.   if (flux_tp == "FluxLinear")
  184.     {
  185.       if (tip.Delta() > grid->Delta(1)/pmin)
  186.     dt = grid->Delta(1)/pmin;
  187.       else
  188.     dt = tip.Delta();
  189.     }
  190.   else if (flux_tp == "FluxSimpleBuckleyLeverett" || flux_tp == "FluxBurger")
  191.     dt = grid->Delta(1) / maxFluxDeriv (flux(),0,uL,0,0,0.1,0.1);
  192.   else if (flux_tp == "FluxBuckleyLeverett")
  193.     {
  194.       real pdelta;
  195.       if (pmin == pmax)
  196.     pdelta = 10;   // void infinite loop in maxFluxDeriv!
  197.       else if (pmax > pmin)
  198.     pdelta = pmax - pmin;
  199.       else
  200.     errorFP("ConsLaw::setIC","pmax=%f < pmin=%f ....",pmax,pmin);
  201.  
  202.       dt = grid->Delta(1) / maxFluxDeriv (flux(),0,uL,pmin,pmax,
  203.                      0.1,pdelta/10.0);
  204.     }
  205.   tip.setTimeStep (dt);    // store the computed dt in the TimePrm object
  206. }
  207.  
  208. void ConsLaw:: timeLoop ()
  209. {
  210.   TRACE("ConsLaw::timeLoop")
  211.  
  212.   tip.initTimeLoop();
  213.   tip.increaseTime();
  214.   while (!tip.finished())
  215.     {
  216.       solveAtThisTimeLevel ();
  217.       saveResults ();
  218.       updateDataStructures ();
  219.  
  220.       tip.increaseTime();
  221.     }
  222. }
  223.  
  224. void ConsLaw:: solveAtThisTimeLevel ()
  225. {
  226.   TRACE("ConsLaw::solveAtThisTimeLevel")
  227.  
  228.   scheme->update ();
  229.  
  230.   // set boundary conditions:
  231.   u->valueIndex(grid->getBase(1)) = uL;
  232. }
  233.  
  234. void ConsLaw:: saveResults ()
  235. {
  236.   TRACE("ConsLaw::saveResults")
  237.  
  238.   // dump plot data values on file casename.pl if any given point of
  239.   // time for plotting is within 0.9*dt from the present time value:
  240.  
  241.   if (time_points_for_plot.within (tip.getTime(), 0.4999*tip.Delta()))
  242.     DrawFD::makeCurvPlot (u(),file,"u(x) for t fixed","u",
  243.               oform("t=%g",tip.getTime()));
  244. }
  245.  
  246. void ConsLaw:: updateDataStructures ()
  247. {
  248.   u_prev() = u();
  249. }
  250.  
  251. void ConsLaw:: resultReport ()
  252. {
  253.   TRACE("ConsLaw::resultReport")
  254.  
  255.   // single line short report of results:
  256.   // area under the u-curve:
  257.  
  258.   int i1 = grid->getBase(1);
  259.   int in = grid->getMaxI(1);
  260.   int i;
  261.   real area = 0.5*u->valueIndex(i1) + 0.5*u->valueIndex(in);
  262.   for (i = 2; i <= in-1; i++)
  263.     area += u->valueIndex(i);
  264.   area *= grid->Delta(1);
  265.  
  266.   s_o << oform("%22s %14s dx=%7.4f, dt=%7.4f, area=%7.4f\n",
  267.            scheme_tp.chars(),flux_tp.chars(),
  268.            grid->Delta(1),tip.Delta(),area);
  269. }
  270.  
  271.   
  272. /* LOG HISTORY of this file:
  273. $Log: $
  274. */
  275.